The concentration of $CO_2$ in the atmosphere has steadily increased since the industrial revolution leading to several serious climate related problems, and consequently an interest in understanding earth's changing climate. Earth's surfaces has the ability to absorb and store $CO_2$ from the atmosphere, playing a crucial role in the carbon cycle on earth. For instance, the global oceans acts as a natural carbon sink that help to buffer the emissions from human activity. This effect has great influence on the climate and it is therefore meaningful to construct scientific models that describes how the ocean reacts to atmospheric changes.
Quantifying the carbon absorbed by the global oceans is a complex problem and to simplify calculations everything is assumed constant in the horizontal direction, which means conditions only depend on depth and time. This is known as the water column model. Furthermore the project model employs the one-dimensional diffusion equation to describe how the concentration of dissolved inorganic carbon (DIC) in the ocean changes. Initially looking at shallow waters and constant atmospheric $CO_2$ concentration and subsequently in deeper water with an annual increase in atmospheric $CO_2$ concentration of $2.3$ppm.
The first sections 2.1 Libraries and 2.2 Constants imports libraries, sets global constants. Values for depth, time and discretization steps for task 1 and task 2 are given in the first part of 3.2 Functions and 4.2 Functions.
Section 2.3 Constructor for L, R and S matrices creates the matrices $L$, $R$ and $S$ as described respectively in equations $(26a)$, $(26b)$ and $(24)$. The function is generalized to work for both task 1 and 2.
The computation of the DIC data for task 1 and 2 is in the sections 3.2 Functions and 4.2 Functions under the headline DIC concentration data. This task is split in two functions where one initiates the matrices and discretization, and calls the other function for the numerical computation. The functions are constructed this way such that the latter function can use the numba library (@njit) for increased efficiency. The time development of the DIC concentration is solved by the matrix equation $(28)$ which is a discretization of the one dimensional diffusion equation $(1)$. The time developed DIC concentration data is stored in a 2-dimensional array named C_data in the functions. Each element of C_data contains the concentration of DIC over depth and corresponds to a time value. Data for DIC, time and depth are stored in task_1 and task_2 for further visualization.
It is important to note that in task 1 and 2 the initial condition and development of $C_{eq}$ and consequently the $S$ matrix is different. In task 1 the zero concentration initial condition has been implemented as a concentration array filled with zeros, and the boundary conditions is given in the $S$ array with a constant $C_{eq}$ as in equation $(24)$. In task 2 the constant initial concentration is given as an array filled with the equilibrium value of $415$ppm. Since $C_{eq}$ varies with each time step in this problem, $S$ is now a matrix filled with the arrays from equation $(24)$ with varying $C_{eq}$.
Furthermore numerical integration of C_data is used to obtain the total mass of DIC and $CO_2$ in the ocean.
Note: $(number)$ represents the corresponding equation from project description [2].
%matplotlib notebook
import numpy as np
from matplotlib import pyplot as plt
from matplotlib.animation import FuncAnimation
from scipy import linalg, sparse, integrate
from numba import njit
from IPython.display import set_matplotlib_formats, display, HTML, Image
# Plot parameters
fontsize = 12
plt.rcParams.update({
"axes.titlesize": fontsize + 5,
"axes.labelsize": fontsize,
"axes.grid": True,
"axes.linewidth": 1,
"lines.linewidth": 1.5,
"lines.markersize": 7,
"figure.figsize": (9, 4.5),
'axes.titlepad': 10,
"ytick.labelsize": fontsize,
"xtick.labelsize": fontsize,
"ytick.major.pad": 5,
"xtick.major.pad": 5,
'legend.loc': 'best',
"legend.fontsize": fontsize,
"legend.handlelength": 1.5,
"legend.frameon": True,
"mathtext.fontset": "stix",
"font.family": "STIXGeneral"
})
set_matplotlib_formats("svg")
# General
k_w = 6.97e-5 # [m s^-1] - Mass transfer coefficient
a_1 = 6.97e-7 # [s m^-1] - Proportionality constant
u = 10 # [m s^-1] - Wind speed
H = 5060 # [mol m^-3 atm^-1] - Proportionality constant
P_CO2 = 415e-6 # [atm] - Partial pressure CO2, assuming ideal gass
C_eq = H * P_CO2 # [mol m^-3] - DIC equilibrium concentration
startyear = 2020 # [year] - First year of simulations
ocean_area = 360e12 # [m^2] - Total area of ocean
C_atomic_mass = 12 # [g mol^-1] - Atomic mass of carbon
# Task 1
K_0_1 = 1e-3 # [m^2 s^-1] - Constant
K_a = 2e-2 # [m^2 s^-1] - Constant
K_b = 5e-2 # [m^2 s^-1] - Constant
z_a = 7 # [m] - Constant
z_b = 10 # [m] - Constant
L = 100 # [m] - Depth of ocean
# Task 2
K_0_2 = 1e-4 # [m^2 s^-1] - Constant
K_1_2 = 1e-2 # [m^2 s^-1] - Constant
a_2 = 0.5 # [m^-1] - Constant
z_0 = 100 # [m] - Constant
def matrix_constructor(
depth, # (int) - Total depth of ocean
N, # (int) - Discretization steps for depth
delta_t, # (float) - Timestep
matrix_name, # (string) - Name of matrix to construct - L, R, S
K_func, # (function) - Function for diffusivity
C_equilibrium = np.array([C_eq]) # (numpy.ndarray) - Equilibrium concentration in ocean
):
### Discretization of diffusivity
z_values, delta_z = np.linspace(0, depth, N + 1, retstep=True)
K_values = K_func(z_values)
### Constants
sign = 1 if matrix_name == 'L' else -1 # L -> 1 | R -> -1
alpha = delta_t / (2 * delta_z**2)
gamma = 2 * alpha * k_w * delta_z * (1 - (K_values[1] - K_values[0]) / (2 * K_values[0]))
### S matrix
if (matrix_name == 'S'):
cont_column_vec = np.insert(np.zeros(N), 0, 2 * gamma)
return np.dot(cont_column_vec[:, None], C_equilibrium[None, :]).T
### L and R matrix
## Main diagonal
main_0 = 1 + sign * (2 * alpha * K_values[0] + gamma)
main = 1 + sign * 2 * alpha * K_values[1:]
main = np.insert(main, 0, main_0)
## Semidiagonals
K_n_marked = K_values[2:] - K_values[:-2]
# Upper diagonal
upper_0 = -2 * alpha * K_values[0]
upper = -alpha / 4 * K_n_marked - alpha * K_values[1:-1]
upper = np.insert(upper, 0, upper_0) * sign
# Lower diagonal
lower_N = -2 * alpha * K_values[-1]
lower = alpha / 4 * K_n_marked - alpha * K_values[1:-1]
lower = np.append(lower, lower_N) * sign
### Total matrix
return sparse.diags([upper, main, lower], offsets = [1, 0, -1]).toarray()
Task 1 looks at the response in shallow areas of the ocean due to $CO_2$ changes in the atmosphere. The problem is of great importance because the increase in DIC concentration leads to acidification which is damaging for wildlife. The Norwegian continental shelf is such an area and the approximated depth is 100m.
The goal is to determine if the changes in the atmosphere and in the ocean are in sync, implying that the DIC concentration is in constant equilibrium with the atmosphere, or if it lags behind.
task_1 = {
'depth': {
'value': 100, # Depth of ocean
'steps': 500 # Discretization steps
},
'days': {
'value': 180, # Days to simulate
'steps': 800 # Discretization steps
}
}
def K_1(
z # (float) - Depth from ocean top
):
return K_0_1 + K_a * (z / z_a) * np.exp(-z / z_a) + K_b * (L - z) / z_b * np.exp(-(L - z) / z_b)
@njit
def get_C_data_1(
z_N, # (int) - Depth discretization steps
time_N, # (int) - Time discretization steps
C, # (numpy.ndarray) - Initial DIC concentration
L, R, S # (numpy.ndarray) - L, R, S matrices
):
C_data = np.empty((time_N, z_N + 1))
C_data[0, :] = C # First element
for n in range(1, time_N):
V = R.dot(C) + S
C = np.linalg.solve(L, V)
C_data[n, :] = C # Insert value
return C_data
def task_1_DIC_data(
depth, # (int) - Total depth of ocean
z_N, # (int) - Discretization steps for depth
days, # (int) - Time in days
time_N # (int) - Discretization steps for time
):
### Discretization
delta_t = days * 24 * 60 * 60 / (time_N - 1)
time_values = np.linspace(0, days, time_N)
z_values = np.linspace(0, depth, z_N + 1)
### Matrices
L = matrix_constructor(depth, z_N, delta_t, 'L', K_1)
R = matrix_constructor(depth, z_N, delta_t, 'R', K_1)
S = matrix_constructor(depth, z_N, delta_t, 'S', K_1)[0]
C = np.zeros(z_N + 1).T # Initial C
### Time development for C
C_data = get_C_data_1(z_N, time_N, C, L, R, S)
return C_data, time_values, z_values
def plot_DIC_C_time_dev_1(
C_data, # (np.ndarray) - DIC concentration data
time_values, # (np.ndarray) - Discretization of time
time_N # (int) - Time discretization steps
):
plt.figure()
plt.title('Time development of DIC concentration')
plt.ylabel('DIC concentration [$mol\cdot m^{-3}$]')
plt.xlabel('Time [$days$]')
# DIC equilibrium
plt.plot(time_values, [C_eq] * time_N, label='$C_{eq}$')
# DIC max
plt.plot(time_values, C_data.max(axis=1), label='$C_{max}$')
# DIC min
plt.plot(time_values, C_data.min(axis=1), label='$C_{min}$')
plt.legend()
plt.show()
def plot_DIC_C_depth_anim(
C_data, # (np.ndarray) - DIC concentration data
z_values, # (np.ndarray) - Discretization of depth
days, # (int) - Time in days
depth, # (int) - Total depth of ocean
anim_steps = 100, # (int) - Animation steps
tracing = False # (boolean) - Trace plot-line
):
plt.ioff() # Turns off interactive plotting
fig, ax = plt.subplots()
plt.title('DIC concentration as a function of depth')
plt.xlabel('DIC concentration [$mol\cdot m^{-3}$]')
plt.ylabel('Depth [$m$]')
fig.set_tight_layout(True)
ax.set_xlim(0, 2.5)
ax.set_ylim(0, 100)
plt.gca().invert_yaxis()
## First line
if tracing:
lines = [ax.plot(C_data[0], z_values, color='b', alpha=.20, linewidth=1, label='Day: 0')[0]]
else:
lines = [ax.plot(C_data[0], z_values, color='b', linewidth=2, label='Day: 0')[0]]
def update2(i):
label = f'Day: {int(np.round(i / len(C_data) * days))}'
if tracing: # See previous lines
label0 = f'Day: 0 - {int(np.round(i / len(C_data) * days))}'
## Label for all previous lines
lines[0].set_label(label0)
## Change previous line style
lines[-1].set_label('_nolegend_')
lines[-1].set_alpha(.20)
lines[-1].set_color('b')
lines[-1].set_linewidth(1)
## Add new line
lines.append(ax.plot(C_data[int(i)], z_values, color='m', linewidth=2, label=label)[0])
else: # Only one line
lines[0].set_data(C_data[int(i)], z_values)
lines[0].set_label(label)
plt.legend(loc='upper right')
return lines, ax
## Create animation
anim = FuncAnimation(
fig, update2,
frames=np.linspace(0, len(C_data) - 1, anim_steps), interval=100,
blit=True, repeat=False
)
## Save animation
return_anim = HTML(anim.to_jshtml())
plt.close()
plt.ion()
return return_anim
data = task_1_DIC_data(task_1['depth']['value'], task_1['depth']['steps'], task_1['days']['value'], task_1['days']['steps'])
task_1['C_data'] = data[0]
task_1['time_values'] = data[1]
task_1['z_values'] = data[2]
plot_DIC_C_time_dev_1(task_1['C_data'], task_1['time_values'], task_1['days']['steps'])
Plot shows extreme values of DIC concentration in the ocean as a function of time. The simulation shows development over 180 days, assumes zero initial concentration of DIC in the ocean and 415 ppm $CO_2$ in the atmosphere.
display(plot_DIC_C_depth_anim(task_1['C_data'], task_1['z_values'], task_1['days']['value'], task_1['depth']['value'],
tracing = True))